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Abstract. Using radioactive nuclei for electricity generation in Microelectromechanical 
Systems (MEMS) is important research as needs for longer life batteries increase. There are 
many applications developed in recent years, however there are limitations still to overcome 
before a final product can be produced. One of the important issue is the low power output. 
This research addresses this issue with a new method in fabrication for powering MEMS sensors. 
We have proposed to fabricate 63 'Ni nano-particle /3~source in a glassy phosphorous type sphere 
which creates scintillation and phosphorescent photons. The micro-spheres will be doped in our 
semiconductor. Since 63 Ni is a pure /3~ emitter, in this report the energy loss dE/dx of f3~ in 
our scintillation material (phosphorus 15 P) is modelled using C++ coding with GEANT4 and 
furthermore the particle distributions in two different source geometries (circular and square 
structure) is studied using Finite Element Analysis (FEA). We have shown the 63iVi j3 energy 
spectrum absorption in 15P, where upto energy range 28keV had nearly 100% efficiency. The 
optimum thickness therefore achieved was 800^im. 



1. Introduction 

Betavolatics are direct charge batteries where charges on (3 particles directly drive a current. In 
this research a new type of Betavolatics are proposed, where a 63 Ni (source) will be fabricated 
inside phosphorus ( 15 -P) spherical containers, and ( 15 -P) elements are then doped in silicon 
lattice. The energy lost from (3 particles to the absorber ( 15 P) is then derived into the solar 
arrays for final conversion essentially the phosphorous-solar cell system acts as an impedance 
matching device converting a high voltage, low current (e.g 17keV single electron) into a useful 
power source (e.g 2V, 100mA. It is convenient to split the energy loss process into two different 
categories on the base of their scattering properties into elastic and inelastic process. In elastic 
scattering the target remains intact and in inelastic scattering the particle loses its kinetic 
energies to the target and causes other internal process. The process of inelastic scattering 
analysed using two different methods, Bethe and Gyrzinski. The Bethe theory is looking at the 
loss by Born approximation to calculate the cross section, where Gyrzinski attempts to obtain 
the cross section by looking at the inner-shells ionisation. The Inelastic Mean Free Path (IMFP) 
is also obtained and reported. The particle flow distributions at two geometries and radiation 
absorption in 15 P are modelled using Neumann-Dirichlet boundary condition. 



2. Theoretical Calculations 

Environmental factors and material surroundings are always denned using boundary conditions. 
Enforcing boundary condition on flux requires an intensive numerical treatment. In this model 
the boundary conditions are defined using Dirichlet equation. By using Dirichlet we are able to 
impose the Laplace equation (our transport equation) to the system domain. Figure [T] demon- 
strates Dirichlet boundary conditions for circle and square geometry for an incoming flux. 




Figure 1. Particles flow on the surface of the 63 source. At figure la enter the surface from 
circle structure where in figure lb enters the surface from a square structure. The unit on the 
source is dimensionless. 



As shown the flux in figure la enters the surface from circle structure where in figure lb enters 
the surface from a square structure. To understand the behaviour first it is more convenient to 
have the numerical Laplace for both cases. 
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Circular Laplace has true Dirichlet boundary everywhere between < r < 1 and U = r m sin{m6) 
PQ . As given in figure [I] the local minima and maxima in the domain is defined on the bound- 
ary (this is one of the Laplace properties as they tend to be harmonic, the value 1.4 is a 
arbitrary threshold that is defined to study the Lagrange function). In the square structure 
we have Dirichlet boundary conditions everywhere except y = 1 (also < x < 1). Thus 
U(x,y) = Y l ^' =1 C n sinh(niry)sin(nirx) and C n = 2 J Q f(x)sin(nnx)dx/sinh(nn). 
Now that Laplace equation is defined we need to numerically define the flux and multiply the 
two values. Thus we will have: 
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Here a is the velocity term, 7 is the source term and equation [2] is called flux vector, c can be 
also indirectly calculated for an anisotropic material. The flux term can be also calculated using 
continuity equation. Hence in more general form (where the material convection, absorption 
and conservative flux convection is set zero), our problem can be simplified and written as: 



(3) 



Equation [3] is in computational domain (uj), thus we need to satisfy all conditions in boundary 
domain. This is called Neumann-Dirichlet where the boundary will be transformed from uj to 
duo (from computational boundary to domain boundary). This transformation is also described 
as domain decomposition preconditioner [2]. Thus the partial differential equation is given by 
V 2 ^ + u = where V donated as Laplacian therefore we will have: 

~T~(. X ) = X7 2 u{x).n(x) (4) 

Where n refers to a normal vector, thus we can rewrite the Eq. 2 as 

n.(c y u + au — 7) = g — h N (5) 



where g and h f N donate the boundary source term and the Lagrange multiplier factor. h% is 
needed in a mixed field situation as it corresponds to local maxima and minima. In some respect 
h l N can also refers to the velocity {3|. Hence we can rewrite the equation [3] as: 

n.( V (p+^)I z )=g-h T N (6) 

This called Neumann generalised boundary condition (in general form the hj^ should be zero. 
Often hjf is replaced by Dirichlet condition and that equals to hu). Now that our Boundary 
Conditions on particle flow are defined, the energy loss properties are also need to be statisti- 
cally tackled. Model I is looking at the Bethe Theory in which the behaviour of each atom is 
studied in terms of transmission from their initial states to final states and the Model II is based 
on the Gryzinski's model looking at the energy loss of atoms due to their inner shell electron 
excitations. As the elastic scattering involves interactions of electrons with atomic nuclei and 
since the atomic nuclei is thousands times more massive than an electron, the energy transfer 
that involves elastic scattering is often insignificant, however this energy can reach up to ten 
thousands of electronvolt for a small fraction of electrons. The general term to determine the 
scattering is given by differential cross section of an incident electron over a solid angle f2: 

da/dn= f 2 (7) 

Where / is the atomic scattering factor which is given in function of scattering angle . The 
differential cross section in elastic form is represented by: 

' \f( q )\ 2 = -±-\Z-f x (q))\ 2 (8) 
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where ao is the Bohr radius 0.529E~ 10 m, Z is the atomic density and the f x (q) is the scattering 
factor for an incident electrons. This usually is set to zero in classical and wave mechanical 
theory. By having the relativistic factor 7 equal to 1 + Eo/(moc 2 ) thus we can have: 

fc - 47222 r i-^™) w 
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where the q = 2kosin(i2) is the momentum transfer; and k is the wave vector of the electron 



before and after the scattering event. 



da 1 = 4fZ 

cM a 2 q± ^ 

Where ro is the screening radius given by Wenz formula [61 in which shows the nuclear potential 
attenuation as function of r distance is given by: 

4> (r) = [Ze/iiTEor] exp (— ) (U) 

The q in form of inelastic will be slightly modified as it depends on the initial and final scatter- 
ing angle. Therefore in inelastic scattering it equals to q 2 = k 2 (6 2 — 6g). Here k = ^ is the 
characteristic angle corresponding to average energy loss. By having the latter in an stationary 
state the % 2 in Eq.ll represents the Rutherford scattering cross section in an elastic form 

a Q q 

[7]. Using Eq.10 and Eq.ll, figure^ is generated to demonstrate the inelastic cross section as 
function of angle. 
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Figure 2. The Angular Dependence of Differential Cross Section of 200eV in 15 P. The inelastic 
scattering cross section decreases as the scattering angle increases 



As Figure [2] demonstrates the inelastic scattering cross section decreases as the scattering angle 
increases. This figure reveals that most of the scattering happens to be in range of Be < < 9q. 
This roughly proportional to I/O 2 . This shows higher probability of scattering in the forward 
direction. Gyrzinski's approach enabled the calculations to measure the energy loss due to the 
excitation. Gryzinskis stopping power simply can be calculated by replacing the [9] calculation 
(Eq.12) by Bethe. Hence the direct simulation of individual excitation specifically in the valance 
band will be more feasible. Since there are a limited number of electrons available in the inner 
shell, the inelastic stopping power is calculated by taking away the portion of energy lost to the 
valance band. Thus we have: 



(12) 



Here the first term is the total energy loss and the second term represent the energy lost to 
valence electrons. Thus i shows the number valance electron which the equation is solved for. 
Now the Gyrzinski Stopping power can be calculated from (for derivation see [9]): 
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Where n% is the number of electrons, Ei is the binding energy and E is the incident energy. 
The second term demonstrates the relativistic factor and the factor 2.7 corresponds to atomic 
ionization potential. Thus the formula for excitation energy is given by: 

d(AE) ~ (AE) 3 E E + Ei [1 E [ E t [1 E ) + ^ J+{ < E% ' 11 

(14) 

AE is the energy loss. Now cross section can be calculated from: 

ai = n ^ 4 i§ ( §rf )L5[1 + 1 {1 ~ i } x ln[2 - 7+ ( l " 1)0 ' 5]] (15) 

Girzinskys formula provides an easy approach to calculate the relative energy loss. The 
simulation is done at 200eV incident electron interacting with three inner shell electrons of 
15 P. Using equation 14 and 15 the magnitude of the excitation energy lost was calculated, and 



that was +69.6eV, hence the due to inelastic scattering is 130. 4eV per 200eV. 
So far the loss due to excitation is determined using Gyrzinski approach thus we required to 
verify the loss by inelastic scattering in more details using Bethe formula. The Bethe theory can 
be used to evaluate the behaviour of each atomic electrons in terms of transition from initial 
state to final state. This process will be angular dependent. flQ] has approached this problem 
using first Born approximation, therefore differential cross section can be determined by: 

da n = (2Tr)~ 2 M 2 fi~' l (k' /k) / exp(iK.r)u n * (n, ...,r z ) x Vu (ri, ...,r z )dri...dr z dr)\ 2 dw, (16) 



where M is the reduced mass of colliding system, r is the position of particle relative to the 
centre of atom, fik is the momentum of particle before the collision and h is the momentum 
after the collision, u is the eigenfunction of rj coordinates whose the total number is Z. Now 
we can have the Eq.2.10 in angular form: 

1^ = ^2^^' / V ( r )^oi'nexp(iq.r)dT\ 2 (17) 



Here the initial states and final states are and ko and K\ are the wave- vectors and 
V(r) is the potential energy required by each collision. The V(r) is required due to existence of 
Coulomb force and it is required by Bethe theory. It can be calculated using: 

Zr 2 1 z p 2 
V (r) = ^--^Y^-^— (18) 
47reo?" Atteq p^\r - Tj\ 

Where the first term represents the Rutherford scattering and |e„(q , )| 2 is the dynamical structure 
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Figure 3. Gyrzinski's differential Cross Section. The differential cross section (in barn) is 
calculated for Gyrzinski (Eq. 14) cross section. Only done for low energy calculation, with 
higher energy considered another rising peak should expected. 



factor. The dynamical structure factor is closely related to the generalised oscillator strength 
|10j . And is calculated using Rydberg energy (13.61ey) and transition energy. Figure [3] demon- 
strates the differential cross section calculated with Gyrzinski's equation (Eq. 14). It is clear that 
the inner shell electrons contribute relatively little to the cross section. One of the major role 
in ionization depends on a process called elastic ionization. To determine the energy loss due 
to elastic scattering the Eq.9 should to be modified. As mentioned before the elastic scatter- 
ing is caused by interaction of incident electron beams with electrostatic field of nuclei. By 
transferring enough energy to the rest mass projectile atom, it will be scattered from the shell 
after the collision and causes the atom to be ionized. In this case the Eq.9 is simply rewritten as: 



da 
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Having electrons scattered elastically makes the second term in Eq.9 zero. This is due to the 
estimation for the screening radius which is defined in Eq.ll. Thus by using Eq.ll one more 
time we can define the differential cross as function scattering angle: 



do 

dn 
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1 is the characteristic angle. The ro can be redefined using Thomas-Fermi 
= ao-Z -1 / 3 ). For an ionized atom the ^ increases with decrease in the char- 



Where #o = (kofo, 
approximation (ro 

acteristic angles. The elastic scattering energy transfer can be also given by: 



E = E max sin 2 {9/2) = E max (l - cos9)/2 = E max (l - cos6)/2 



(21) 



Where E max is the energy transfer at maximum angle and is calculated from [TT] as: 

E 2 max = 2E (E + 2M c 2 ) 



(22) 



The differential cross section is at its minimum at central impact 6 = and increases by increase 
in the angle. 

Now that the loss due to ionization and excitation is measure the IMFP can be calculated. The 
model describing IMFP is based on a theoretical approach to the transport of primary electron 
beams on the surface, thus a relationship between signal intensity and concentration of given 
elements in an elastic and inelastic deferential cross section is required. This could be a complex 
function as the accurate data from elastic and inelastic cross section is not obtainable. Hence 
the computational functions are usually solved by looking into continuous slowing down approx- 
imation (CSDA) in which an electron energy along the trajectory is assumed to be a function 
of length. To tackle this problem the excitation function for a charged particle is defined by its 
dielectric function. As demonstrated by [12J: 
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d(AE)dq = ^tfNE { ^q) ] ~q 1 ' 

Where N is number density of the atom and E is the energy of incident beam, (AE) is the 
energy loss and Jm( £( .~ 1 (? . ) ) is the dielectric function. Dielectric function can be reduced to a 
more generalised form of Lindhard type dielectric function [13] and [H]. In This model the 
IMFP is calculated based on [15] and [H] approach: 

E 

A ~ E 2 [pin(0. 191 p-°- 5 E) - (1.97 - 0.91(N v p/M)/E) + ((53.4 - 20.8{N v p/M)/E 2 )} U ' ' 

Here the p is the material density and E g is the band gap energy, M is the atomic weight, 
f3 = -0.1 + 0.9M(E 2 + E 2 )~ - 5 + 0.069p 01 ), E P = 28.8(N v p/M) - 5 is the free electron plasmon 
energy, and N v is the number of valence electron per atom. 
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Figure 4. Inelastic Mean Free Path in 15 P. As the incident energy increases the IMFP also 
increases. 



Figure |4| demonstrates use of Eq. 2.19 for 15 P in our case. In this case the N v is set to 5, 
M = 30.97, E q = 3eV and p = 1.88g/cm 3 . These are the value which are typical for our 
material (data from NIST library). 




Figure 5. The Particle Distribution in Ni. The distribution in circular geometry is lower in 
centre. However the particles are equally distributed on the surface whereas the square shape 
has a higher concentration in centre and more around the edges with less equal distribution on 
the surface 



3. Results 

Subatomic particle flow distribution is an important factor in designing and manufacturing of 
nano devices. This is vital where source geometry can make a significant difference in making 
good use of particle flow. In this research two geometry system is considered square and circular 
shape. Figure [5] demonstrates the two distributions for a circular and square shape source, the 
flow in circular geometry is more distributed on the surface whereas the flow in square structure 
is piled at the corners. Likewise the edges around the square shape caused more concentration 
in the centre. By taking away the flow in circular form from square shape figure [6] has achieved. 
The geometric flow properties of particles, effects their behaviour due to their interactions and 
their collision dynamics. This could be exhibited by drag force (this is equivalent for electrons), 
the boundary layers on the surface and wake vortices. Also by long and short range inter-particle 
forces. Our model tends to have a Normal distribution (Gaussian Distribution) in space. As 
shown in figure[5]the distribution is balanced in circular shape, where as in square form, particles 
are mounted around the corners. The Gaussian shape has following property: 



1.4e+04 
1.2e+04- 
le+04- 

(fi 

1 8,000 
o 

o 

6.000 
4.000 
2,000 

o - 



Green: Circular 
Blue: Square 
— Gauss Fit 




200 300 
Width 



Figure 6. The Particle Distribution in 63 iVi. The distribution in circular geometry is lower in 
centre. However the particles are equally distributed on the surface whereas the square shape 
has a higher concentration in centre and more around the edges with less equal. The sharp 
peaks at the centre and the corners are also representing noise in the system. 
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Here W is the width, x is the length, A is the area and yo is the offset. The x 2 /DoF calculated 
was 9.13(+5). The black area is the area overlapped for both circular and square shape source. 
As the absorption is concerned, the main processes involves in this model are excitation, ion- 
ization. As figure [7] denotes, at energies over 28keV the incident particles start to deposit most 
of their energies deeper in the film (800 [im thick) with some also penetrating, this means the 
IMFP is increases. This also corresponds to higher particle range. 

The energy loss in 2keV to 28keV range could be in form of heat or bremsstrahlung where the 
loss is negligible and they deposit almost 99.9% of their initial charge, and 1 % is where the 
particles break out the film in form radiation. From 28keV to 60keV range a portion of incoming 
particles escape from the film. However there will be some losses but particles tend to penetrate 
more into the film before their first collision occurs. Thus Figure 7 presents the fruit of our 



calculation and model, it shows the 63Ni f3 energy spectrum absorption in 15 , and the linear 
line corresponding to full energy deposition (100% efficiency) 
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Figure 7. The Energy Absorbed in 15 P. The film of 15 P is set to have 800 \im. A monoenergetic 
particle source is created with different energies between E m i n = 2keV and E max = QOkeV to 
simulate the 63 Ni spectrum. As figure denotes, at energies over 28keV the incident particles 
starts to deposit most of their energies deeper in the film with some also penetrating, this means 
the IMFP is increases. This also corresponds to higher particle range 

Particle range consists of four profiles, longitudinal range, lateral straggling, projected range, 
and collisional range. The longitudinal range increases because of increase in energy spread(of 
incoming beam), which is not influenced by dispersion effect. This has a direct relation with 
beam profile. Although this could be sometimes negligible particularity for a single beam shot 
but cannot be ignored for a continuous beam (primarily for lower energies). Considering the 
initial beam profile be e z the final Root Mean Square (RMS) is given by [17J: 

elfinal =4 + (^) 2 (X' 2 ) 1(Z 2 ) + °?D 2 (5 2 ) + 2aD (Z5)} (26) 

Where aD is the longitude dispersion, D is the dispersion factor and X,Z,5 are the angle, 
position and energy spread respectively. More studies can be done during the experimental pro- 
totyping and the effect of beam dispersion could be taken into account before a valid comparison 
between the final prototype and simulation can be made. 

Another dependent of particle range is lateral scattering. As demonstrated in figure [8] beam 
profile does not significantly change if the distance to any field edges is greater one half of elec- 
tron ranges [18] . Lateral straggling is significantly affected by size of incoming beam and this is 
because of correlations between the particle range and its position. 

The projected mean range is calculated from CSDA function. The simplified mean free particle 
range hence is given by [19J: 



(R) = t(l + t/X ) 



(27) 



Figure 8. Particle Range Profiles in Phosphorus. The green line demonstrates the Lateral 
range. Lateral range is largely effected by the beam size, and this is because of particle 
correlations between range and their positions. The red line shows the longitudinal range. 
The longitudinal range occurs from beam energy spreading (which is not effected by dispersion) . 
The blue line demonstrates the projected range and can be given using CSDA and finally the 
red-star line shows the stopping range where the particles deposit most of their energies in the 
bulk. The red line is the one we are most interested in as a first order approximation 



Where Ao is mean free path of electrons, and limited to value of thickness (t). When t is ap- 
proaching R (range), the corresponding value of R can be called L (mean free range), thus L 
can be recalled as: 

fl = A [(4L/A) + 1 f-l 



Here A = Ao/4 and can be given as [20] : 

X = 2a 1 (Z)E^ {z) (29) 

In experimental approximation a n {z) is the tabulated thickness where n = 1,2, .... Finally the 
collisional range is where the incoming beam deposits most of their energy due to the collision. 
This value represents an accumulated depth in the film where most of the interactions take 
place. This could later facilitate an understanding of an optimum thickness required for the 
application. 



4. Conclusion 

This report has discussed modelling of 63 iVi /3-source doped in phosphorous ( 15 -P). The 
energy loss processes are discussed and analysed theoretically and experimentally. The Q value 
calculated for our (3~ source ( 63 iVi) was 67keV with mean energy of 17keV. Interactions of 
f3 particles with 15 P bulk has led to energy loss and charge deposition. Thus to measure the 
maximum output charge first it is appropriate to look at the number of decays per second. Hence 



having 63 Ni with half life of 3.1567(+9)s, the probability of decay in one second is 2.196(— 10). 
Also by having mass of 63 Ni as 1.044(— 22)g the number of atoms in 10mg will be 9.57(+19) 
hence n is equal to 2. 103 (+11) (this values represent the theory, of-course the real application 
values are much less). The energy loss analysed in this model by implementing discrete points 
from 63 Ni — f3~ spectrum to simulate the interaction of incoming particles. The total energy 
loss calculated in 800//m 15 P was 5.63(+5)eV.s _1 . However a portion of incoming beams will 
be lost due to penetration through the film. 
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